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The stability and chaotic behaviors of Bose-Einstein condensates with two- and three-atom in- 
teractions in optical lattices are discussed with analytical and numerical methods. It is found that 
the steady-state relative population appears tuning-fork bifurcation when the system parameters 
are changed to certain critical values. In particular, the existence of three-body interaction not only 
transforms the bifurcation point of the system but also affects greatly on the macroscopic quantum 
self-trapping behaviors of the system associated with the critically stable steady-state solution. In 
addition, we also investigated the influence of the initial conditions, three-body interaction and the 
energy bias on the macroscopic quantum self-trapping. Finally, by applying the periodic modulation 
on the energy bias, we find that the relative population oscillation exhibits a process from order to 
chaos, via a series of period-doubling bifurcations. 
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I. INTRODUCTION 

In Recent years, Bose-Einstein condensates (BECs) in optical lattices have attracted enormous attention both 
experimentally and theoretically [J, Q • This is mainly because the lattice parameters and interaction strength can 
be manipulated using a modern experimental technique. Making use of this, researchers have discovered many long- 
predicted phenomena, for example non-linear Landau-Zener tunneling, energetic and dynamical instability and the 
strongly inhibited transport of one-dimensional BEC in optical lattices 0,11, IE B SJE H More attracting 

phenomena, namely, self-trapping, was recently observed experimentally in this system [111 ]. In such an experiment, 
a BEC cloud with repulsive interaction initially loaded in optical lattices was self-trapped. Many theoretical analysis 
was also presented about self-trapping [l2|, [HI, EH, [H| • It is well know the macroscopic quantum self-trapping (MQST) 
means self-maintained population imbalance with non-zero average value of the fractional population imbalance which 
was detailed discussed [lg, [TtJ • Marino et. al. considered that the damping decays all different oscillations to the 
zero-phase mode [Hj]. Besides, macroscopic quantum fluctuations have also been discussed by taking advantage of 
second-quantization approaches fl9|. However, when the trapping potential is time dependent and the damping 
and finite-temperature effect can not be neglected, chaos emerges. Abdullaev and Kraenkel studied the nonlinear 
resonances and chaotic oscillation of the fractional imbalance between two coupled BEC's in a double- well trap with 
a time-dependent tunneling amplitude for different damping [20]. When the asymmetry of the trap potential is 
time-dependent and its amplitude is so small that can be took as a perturbation, Lee et al. studied the chaotic and 
frequency-locked atomic population oscillation between two coupled BECs with a weak damping, and discovered that 
the system comes to an stationary frequency- locked atomic population oscillations from transient chaos (2lj . 

It is important to note that theoretical studies of stability are mainly focused on the effect of two-body interactions. 
It is clear that in low temperature and density, where interatomic distance is much greater than the distance scale of 
atom-atom interactions, two-body s-wave scattering should be important and three-body interactions can be neglected. 
But, if the atom density is higher, for example, in the case of BEC in optical lattices, three-body interactions will 
play an important role |22] |. As reported in Ref. (23|, even for a small strength of the three-body force, the region of 
stability for the condensate can be extended considerably. 

Therefore, the purpose of this paper is to investigate the steady-state solution of BEC in an one-dimensional 
periodic optical lattice when both the two-body and three-body interactions are taken into account. By using the 
mean-field approximation and linear stability theorem, one interesting result is found that the tuning-fork bifurcation 
of steady-state relative population appears when the system parameters are changed to certain critical values. The 



* Corresponding author. Email: ychen@gmail.com 



2 



existence of three-body interaction not only transforms the bifurcation point of the system but also affects greatly 
on the self-trapping behaviors of the system associated with the critically stable steady-state solutions. Additionally, 
we also study the effects of the initial conditions, three-body interaction and the energy bias on the MQST. Besides, 
we discuss the chaos behaviors of the system by applying the periodic modulation on the energy bias. The result 
shows the relative population oscillation can undergo a process from order to chaos, via a series of period-doubling 
bifurcations. 

This paper is organized as follows. In Sec. II, we introduce the mean-field description of BEC in optical lattices with 
two- and three-atom interactions. In Sec. Ill, with linear stability theorem, we analysis the stability of steady-state 
solutions. Then the influences of three-body interaction on the macroscopic quantum self-trapping of the system are 
displayed In Sec. IV. In Sec. V, by applying the periodic modulation in the energy bias, we discuss chaotic behaviors 
of the system using the numerical simulation method. In the last section, summary and conclusion of our work are 
presented. 



II. MEAN-FIELD DESCRIPTION OF BEC IN OPTICAL LATTICES WITH TWO- AND 

THREE-ATOM INTERACTIONS 

We focus our attention on a BEC with both two- and three-body interactions is subjected to one dimensional 
(ID) optical lattices where the motion in the perpendicular directions is confined. In the mean-field approximation 
, the dynamics of BEC can be modeled by the ID Gross-Pitaevskii (GP) equation in the comoving frame of the 
lattice 0,0,0 HI, 

lh i^ = ~^L{ h §i~ imait ) ^+ v oco S (2K lX )^ + ^^m^ + ^- r m^, (i) 

where $ is the wave function of the condensate, m is the mass of atoms, a s is the two-body s-wave scattering length, v o 
is the strength of the periodic potential, Ki is the wave number of the laser light which is used to generate the optical 
lattice, mai stands for either the inertial force in the comoving frame of an accelerating lattice or the gravity force, 
a± = y 1 'h I '(mui ±) , where ui± is the radial frequencies of the anisotropic harmonic trap, (72 |<I?| 4< J?/ (37r 2 a^_) is three-body 
interactions related to the GP equation. Among Eq. ([T]), all the variables can be rescaled to be dimensionless by the 
following system's basic parameter x ~ 2Kix, $ ~ ^/ag JV ' * ~ m^f^' we obtain the normalized 1D-GP equation in 
optical lattices with cubic and quintic nonlinearities, 

i^- = -1 - iat\ $ + ucos(a;)$ + c|$| 2 $ + A|<I>| 4 $, (2) 

where v — ^i^-i , a — 8 ^ K n ai , c = ^~{t2 is the effective two-body interaction, N is the total numbers of atoms, 

A — 3 ™2^2^4 is the effective interaction among three atoms, here the three-body interaction is expected to be positive 
with a value of < A < 1. 

In the neighborhood of the Brillouin Zone edge k — 1/2, the wave function can be approximated by [H 

${x, t) = a{t)e ikx + b{ty {k ~ l ^ x , (3) 

where a(t), b(t) are the probability amplitudes of atoms in each of the two wells respectively and \a\ 2 + \b\ 2 = 1. By 
inserting such wave functions into Eq. ^ and performing some spatial integrals, we obtain the dynamical equations 
with two- and three-body interactions. 

= la + \ (|6| 2 - H 2 ) a + A (l + 2|a| 2 |6| 2 + 2|6| 2 ) a + V -b, (4) 

4; = -^-|(H 2 -H 2 )6 + A(l + 2H 2 H 2 + 2H 2 )6+^a. (5) 

Here, the level bias j(t) = at, and a is the sweeping rate, c and A represent the nonlinear parameters, v is the coupling 
constant between the two condensates. We introduce the relative population variance 

S = \b\ 2 -\af, (6) 

with the parameters a = \a\e l 9a, b = \b\e l 9b, 

9 = 9 b -9 a . (7) 
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Combining Eqs. (|4][7]), one yields the equations of the relative population and relative phase, 



s = —vy 1 — s 2 sin#, (8) 

US 

9 = 7+(c + 2A)s+ cos6>. (9) 

VI — s 2 

s and 9 denote the time derivative of the relative population and the relative phase. If we regard s and 9 as the 
canonically conjugate variables Eqs. (8) and (9), become a pair of Hamilton's canonical equations with the conserved 
effective Hamiltonian 

H = 7S+ i(c+ 2A)s 2 + v\fl - s 2 cos 9. (10) 

In the following section, we will discuss the stability of steady-state in the symmetric condition (7 = 0) with linear 
stability theorem. 

III. STABILITY ANALYSIS OF THE STEADY-STATE SOLUTIONS 

In Sec. II, we have given the dynamical equations of the system with three-body interaction. In this section, we will 
discuss the stability of steady-state in the symmetric condition. Generally, there are two ways to study the stability 
of nonlinear system, the linear stability theorem and the Lyapunov direct method. We will investigate the stability 
of the system with the first method. 

The steady-state solution of this system can be obtained by setting Eqs. (8) and (9) to zero. The forms of steady- 
state solutions are very complicated when the level bias 7 7^ 0. For simplicity, we set 7 = 0, leading to 



s 



h{s,9) = -v^/l-s 2 sine, (11) 

(12) 

(13) 



= f 2 (s,9) = {c + 2\)s+ VS cosO. (12) 

V 1 — s 



and the conserved energy 



Taking s = 0, = 0, we get 



H = i(c + 2A)s 2 +vy/l- s 2 cosO. (14) 



-v\J\-s 2 sin0 = 0, (15) 

US 

(r.+ 9\)»+ rnsfl = 0. (16) 

yl — s 2 



The steady-state solutions obeyed Eqs. (14) and (15) regard as 



9i = 2mr, si = for H = -v, (17) 
h = (2n+l)7r, s 2 = for H = v, (18) 



^—) 2 for H=^l±f 
y c + 2\' 2(c+2\y 



73.4 = (2n + 1)tt, S3i4 = ± a /1-(_- t )2 for H = v o/ _ / n . (19) 



According to the linear stability theorem, we look for the perturbed solutions which are near the steady-state solutions, 

fl(t)=«<(t)+ei(t), 6(t) =6i(t)+c 2 {t) (20) 

where Si(t), (£) for i = 1,2,3,4 signify the steady-state solutions, |ei(i)| -€i \si(t)\ and | £2 ) | <§C \ei{t)\ which is 
relate to the first-order perturbed. Inserting the above expression into Eqs. (11) and (12), we can obtain the linear 
equations near to the steady-states of the nonlinear equations as 

6l= (~is^) £l + (l^O £2 namel y e\ = anE! + a 12 £2 (21) 
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ds J 2 \ dO 7 2 



e'a = ( ~Jr- ) ei + ( "T^r ) £2 namely £2 = a 2 i£i + a 2 2£2 (22) 



such that 



Now, we make use of the above expression to investigate the stability of the steady-states of Eqs. (16-18). 

(l)For 9\ — 2mt,s\ — 0, H — —v, we can calculate the matrix elements an = 0, ayi — —v, 021 = (c+ 2A) + u, 

— u 

a 22 = 0. So, the coefficient matrix of the linearized equations (20) and (21) becomes Ai 



the characteristic equation writes det(Ai — XI) = 



c + 2\ + v 
= 0, which reveals that A 2 + v(c + 2A + v) = 0. 



- A -v 
c + 2X + v 0-A 

We solve the equation to get the two eigenvalues of the matrix A as Ai = y/—v(c + 2A + v), A 2 = —y/—v(c + 2A + v). 
In response to the forms of the eigenvalues, there exist two cases for the stabilities: 

(a) v(c + 2A + v) > 0, that is 

v > and (c + 2A) > -u (23) 

v < and (c + 2A) < -u (24) 

so the two eigenvalues are both pure imaginary numbers. Thus, the stability of the steady-state solutions (0i,si) 
corresponds to a critical case [26j and the dynamical bifurcations between the unstable and stable steady-states will 
appear when the parameters with two- and three-body interactions are changed. 

(b) v(c + 2 A + v) < 0, namely 

v > and (c + 2A) < -v (25) 

v < and (c + 2A) > -v (26) 

so the two eigenvalues are real number. It means that E\ and £2 tend to infinity with the increase of time, and the 
steady-state solutions (#1, s\) are unstable. 

(2)For 82 = (2n+l)7r, S2 — 0, H = v, the matrix elements write as an = 0, 012 = — f, 021 = (c+2A)— v, 022 = 0. The 
corresponding eigenvalues of the matrix A2 become Ai = y/—v(v — (c + 2A)), A2 = —y/—v(v — (c + 2A)). Similarly, 
there are two cases of the stabilities: 



(a) v(v - (c + 2A)) > 0, that is 



(c + 2A) > and v > (c + 2A) (27) 
(c + 2A) < and v > 0. (28) 



so the two eigenvalues are both pure imaginary numbers. And the stability of the steady-state solutions (62,82) of 
the nonlinear equations are reviewed as critical and the dynamical bifurcations will occur, 
(b) v(v - (c + 2A)) < 0, that is 

v > and (c + 2X) >v (29) 
(c + 2A) < v and v < (30) 

so the two eigenvalues are positive or negative real number, respectively. £1, £2 tend to infinity as increasing the time 
to infinity, and the steady-state solutions (8 2 , s 2 ) are losing their stability. 

(3)For 8 3A = (2n+l)7r, s 3i4 = ±Jl — ( c +2a ) 2 ' H ~ ^2(0+2 a)" ' tne ma t r i x elements read an = 0, a i2 = v 2 / (c+2A), 



021 = (c + 2A) - (c + 2A) 3 /u 2 , a 22 = 0, and the eigenvalues Ai = ^v 2 - (c + 2A) 2 ), A 2 = -\Jv 2 - (c + 2A) 2 ). In 
Eq. (18) the population 83,4 are both real quantities which implies 

(c + 2A) 2 >u 2 (31) 

Therefore, the two eigenvalues are pure imaginary numbers. The stability of the steady-state solutions (#3,4, S3. 4) of 
the nonlinear equations are regarded as critical and the dynamical bifurcations will emerge at the bifurcation point 
(c + 2A) = v, s = 0. Obviously, the existence of three-body interaction can change the bifurcation point of the system. 
It plays a important role for stability analysis of the system, as shown in Fig. Q] For C ™ 2X > 1, the system is in 
the critically stable steady-state (9%, 8%), and for < 1> ($2,s 2 ) is unstable and the two steady-state solutions 

(#3,47 53,4) are critically stable. This is a typical tuning-fork bifurcation, and the bifurcation point is c + 2 \ = 1 

According to the above analysis, we conclude that three steady-state solutions possess different stability for different 
parameter regions. And it is very interesting to arrive at the critically stable steady-state solution in experiment which 
relate to the stable stationary MQST [26]. In the following section, we will illustrate the MQST of the non-stationary 
states in detail by two different methods. 
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FIG. 1: Plots of the tuning-fork bifurcation from Eqs. (17) and (18), where S2,S3,S4 are the steady-state solutions and the 
bifurcation point is C ^ 2X = 1 



IV. THE MACROSCOPIC QUANTUM SELF-TRAPPING OF BEC WITH TWO- AND THREE-ATOM 

INTERACTIONS 



In this section, we investigate the macroscopic quantum self-trapping by plotting the phase trajectories and the 
time evolution of the relative population of the system. 



A. The phase trajectories diagram 

The macroscopic quantum self-trapping refers to the phase space trajectories whose the relative population is not 
equal zero. This can be well understood from the analysis Eqs. (8)-(10), corresponding to the critically stable steady- 
state solutions discussed in sec. II. Three kinds of cases occur with different three-body interaction parameters, as 
shown in Fig. 2. 

(1) In the case of v = 0.2, c = 0.1, < A < 0.05 in the phase space , there are two stable points Pi, P% at s = 0, 6 = 7T 
and s = 0, 9 = respectively [Fig. 2(a)], from the circumstance described by Eqs. (22) and (26). Obviously, for the 
stable points Pi, P2, the atoms distributions are equal in the two adjacent wells, the relative population of the 
trajectories around them is equal to 0. It means that atoms oscillate between two adjacent wells and the macroscopic 
quantum self-trapping phenomenon does not emerge in this case. 

(2) When parameter is set to v — 0.2, c = 0.1, 0.05 < A < 0.15, two more fixed points emerge in the line 9 = n 
marked by P3, P4. Among them, Pi, P3 are steady which is corresponding to condition of Eq. (30). They are located 

in ,s = ±.y/l — ( c + 2 \ ) 2 i uence i P4 is unstable point which lies in s = and corresponds to condition of Eq.(26). As 
seen from Fig. 2(b), for the stable points Pi,P3, the atoms distributions are not equilibrium between two adjacent 
wells, and the relative population of the trajectories around them is not equal to 0. It indicates that atoms are 
self-trapped in one well. We take it as oscillating-phase-type because the relative population s and the relative phase 
9 oscillate around the fixed points. 

(3) For v — 0.2, c = 0.1, A > 0.15 , It emerges new trajectories , i.e. the trajectories across point P c [Fig. 2(c)]. Only 
the fixed point P2 is stable which is relate to Eq. (22). So for these trajectories, s varies with time from region of [— 1, 0] 
to [0,1], Apparently (s) 7^ 0, atoms are self-trapped in one well. We regard it as running-phase-type macroscopic 
quantum self-trapping, as described in Refs. [27l. |28|] and observed in experiment .29]. 

The above changes on the topological structure of the phase space are concerned with the change of the energy 
profile. When the relative phase is zero or tt, energy relying on the parameter with three-body interaction and the 
average population s can be derived from Eq. (10). Seeing Fig. 2 , the transition from case(l)to case(2) corresponds 
to the bifurcation of the energy profile of 9 = tt: energy curve bifurcates from a single minimum to the curve of two 
minima. It means the system goes from the Rabi regime into the self-trapping regime through this bifurcation. The 
lowest order of energy profile with 9 = is — c+ £ x , and the energy of the unstable point P4 is — v which is located on 
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FIG. 2: Trajectories on the phase space of the system with three-body interaction varying from A = to A = 0.25(the first 
row) . Corresponding to in the second row we plot the energy profiles for the relative phase 8 = (red dashed) and 6 — n (blue 
solid) 



the maximal order of energy profile with 8 = it. The results displayed by the phase space trajectories conform to the 
case of steady-state solutions discussed in Sec. III. The transition from case (2) to case(3) is signified by the overlap 
of the two energy regions of the profile. In this condition the trajectory stared from s = — 1, 9 = should be confined 
to the lower half of phase plane, corresponding to the running-phase-type macroscopic quantum self-trapping. 

Connecting the analysis of the steady-state solutions to the above analysis on the energy profile, it concludes that 
stable behaviors of the system change constantly with the increase of A and we obtain a general criterion for the 
macroscopic quantum self-trapping trajectories, namely, H (s, 9) < —v. It plays a critical role to find the transition 
parameters of macroscopic quantum self-tapping. 



B. Numerical simulations of the MQST 



Now, we focus on the dynamic behavior which dominated by Eq. (8) and (9) without the time-dependent system 
parameters. We study the effect parameters of the system on the MQST with numerical method starting form Eq. (8) 
and (9). 

Choosing initial condition s(0) = 0, 9(0) — ir/2, the time evolutions of the relative population Fig. (3a)-(3d) show 
some very absorbing features. In Fig. 3(a), the oscillations are regular and the average the relative population s is 
zero for symmetric well case (7 = 0) with a special parameter, but the corresponding MQST does not appear. If we 
increase A from 0.45 to 0.95 in Fig. 3(b), the MQST does not still appear, but the oscillating period becomes short. 
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FIG. 3: The time evolution of the relative population from Eqs.(8) and (9) with initial conditions s(0) = 0, 8o = n/2 and 
parameter: (a) c = 0.1, A = 0.45, v = 0.2, and 7 = 0; (b) c = 0.1, A = 0.95, v = 0.2, and 7 = 0; (c) c = 0.1, A = 0.45, v = 0.8, 
and 7 = 0; (d) c = 0.1, A = 0.45, v = 0.2, and 7 = 0.5; (e) c = 0.1, A = 0.95, v = 0.2, and 7 = 0.5; (f) c = 0.1, A = 0.45, 
v = 0.8, and 7 = 0.5; 



Similarly, rising v , we obtain the same result as shown in Fig. 3(c). 

Here, we study impacting asymmetric well case (7 ^ 0) on the MQST. when we enhance the level bias to 7 = 0.5 the 
average the relative population is changed to —0.41 in Fig. 3(d). Correspondingly, the oscillating period of s is longer 
and the MQST emerges. Note that parameter c, A and v impact greatly on the MQST which arc plotted in Fig. 3(e) 
and (f). In fig. 3(e), when A is from 0.45 to 0.95, the MQST is suppressed with shorter oscillating period. Similarly, 
with increasing v, the average relative population are changed to —0.21 and the oscillating period becomes shorter 
again, as seen in Fig. 3(f). Thus, the influence of parameter c ,A ,v and 7 on the MQST of the system is very dramatic. 
In the case of 7 = 0, fixing the other parameters and changing the initial condition from s(0) = 0, 0(0) — n/2 of Fig. 3 
to s(0) = 0.8,0(0) = § and s(0) = 0.8,0(0) = tt, we observe that the MQST always emerges with varying s(O),0(O). 
The oscillating period is decreased comparing to Fig.3(a)and Fig. 3(d), but the s is increased to —0.86, —0.72 as shown 
in Fig. 4. 

According to the above analysis, we can draw conclusion that when the initial conditions s(0) = 0, 0(0) = tt/2 are 
read, the parameter c, A, v can impact on the MQST for asymmetric well case(7 7^ 0). In addition, in the symmetric 
case, the MQST does not appear and those parameters only affect the oscillating period of the system. Besides, the 
initial conditions can impact the MQST for anyone parameter set. 



V. NUMERICAL SIMULATION OF CHAOS BY APPLYING PERIODIC MODULATION ON THE 

LEVER BIAS 

As a whole, the elementary features of chaos is that the dynamic behaviors are unpredictable for a deterministic 
system. It is very sensitive for the initial conditions and parameters of the system. So, according to these characteris- 
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tics, we can adjust the parameters to make the system get into or get out of the chaos, in other words, we can control 
the regime appearing chaos. In this section We discuss the chaotic behaviors of the system by numerical method. 

If we apply periodic modulation on the lever bias 7 = Aq + Axsin(u>t), the chaos will appear in a special region, 
where Ao,Ai stand for initial phase and amplitude respectively. Inserting this into Eqs. (8)and (9), one derives the 
below dynamic equation. 

s = -vy/l - s 2 sin6> (32) 

US 

9 = A +Aisin(ujt) + (c + 2\)s+ -==005 (33) 

Vl - s 2 

Starting from Eqs. (32), It is found that the dynamics behavior of the system is periodic in some special parameters 
region and it will vary from order to chaos with the increase of A\ , as shown in Fig. 5. With initial conditions 
s(0) = 0, 9(0) = 7T, the phase orbit is a period-one cycle and the corresponding oscillation is a Rabi oscillation for the 
set of parameters with amplitude Ai — 0.002, as in Fig. 5(a). In this case, we set the oscillating period of the relative 
population T. When A\ = 0.009, the phase orbit becomes period-two in Fig. 5(b). It means the oscillating period of 
s arriving at 2T. Then the phase orbit increases from that of period-four to period-eight with increasing A\ as shown 
in Fig. 5(c)and (d). Fig. 5(e) and 5(f) are plotted for A\ — 0.3 and Ai = 1, where the phase orbit does not show a 
clear periodicity which signifies the emergence of chaos. 

From the above analysis, we find that the oscillating period of the relative population varies from a period-one 
limit-cycle to period-two to period-four and then to period-eight and finally all limit-cycles tend to infinity with 7 
increasing. It exhibits a process from order to chaos, through the period-doubling bifurcations [26]. That is to say, for 
a set of fixed parameter v, c, A, Aq, Ai, s(0), 8(0) and uj, the first-order derivative of relative population transform 
from the single period to multiple period and get into chaos at last with the increase of vibration amplitude A\ . 

For the aim of showing the chaotic MQST, we present the plots of the time evolution of the relative population 
and corresponding plots of power spectra from Eqs. (31) and (32) in Fig. 6. And the parameter of Fig. 5(a) is accord 
with Fig. 6(a) and 6(c) where the system oscillates periodically. Making use of those parameters of Fig. 5(e), we 
plot Fig. 6(b) and 6(d). It shows that the power spectrum appears confusion and the average value of the relative 
population is less than zero, which implies the existence of the chaotic behaviors . 



VI. SUMMARY AND CONCLUSION 



In this paper, we study the stability and chaos of BEC with repulsive two- and three-body interactions immersed 
in a one-dimensional optical lattice. The stability of the steady-state solution are analyzed with the linear stability 
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FIG. 5: Dynamical phase orbits of the dimensionless variables (s,ds/dt) from Eqs. (31) and (32) with parameters v = 0.001, 
c = 0.1, A = 0.45, to = 0.1, s(0) = 0, 6 = tt, and (a) A 1 = 0.002, (b) A l = 0.009, (c) A Y = 0.04, (d) A 1 = 0.12, (e) Ai = 0.3, 
(f)=Ai = 1. Here, Ai denotes the amplitude of the time-dependent relative energy. 



theorem. The analytical results show: (1) For v > and c + 2A > — v or v < and c + 2A < —v, the stability of 
the steady-state solution^ = 2nir, si = 0) is in the critical case. (2) For c + 2A > and v > c + 2A or c + 2A < 
and v > 0, the steady-state solution(#2 = (2n + 1)7T, S2 = 0) is the critical stability. (3) For (c + 2A)2 > t>2, the 
steady-state solution (6^4 = (2n + 1)7T, £3,4 = ±y^l — ( ^a ) 2 ) ^ s a ^ so critically stable. When these relationship are 
not satisfied, the corresponding steady-state solution are unstable. A typical tuning-fork bifurcation of steady-state 
relative population appears in special parameter region. And the existence of three-body interaction can change the 
bifurcation point of the system, which is shown as Fig. 1. It plays a important role for stability analysis of the system. 

The critically stable steady-state solution indicates the existence of the stationary MSQT. The stable behaviors 
of the system change constantly with the increase of A and get a general criteria for the self-trapping trajectories, 
H < —v. In addition, we also investigate the effects of the initial conditions, a set of parameters c, v, A, 7 on MQST. 
It shows that c, v,X could affect on the MQST when s(0) = 0, &o — tt for 7 7^ 0. Particularly, the initial value 
s(0) = 0, 6*o = tt or s(0) = 0, 80 = tt/2 can directly impact on the MQST. Finally, we discuss the chaos behaviors by 
applying the modulation on the energy bias (7 = Aq + A\sinujt). In this case, the system will go into chaos through 
the period-doubling bifurcations with the increasing of A, and the time evolution of the relative population and power 
spectra indicate the existence of the chaos MQST. It suggests that one can adjust the lasing detuning and intensity to 
change the values of the parameters in experiments. This adjustable parameters supply the possibility for controlling 
the instabilities of the system, MQST state and the chaotic behaviors. 
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Frequency(z) Frequency (Hz) 



FIG. 6: (a) and (b): The time evolution of the relative population of the relative population from Eqs. (31)and (32) with the 
parameters v = 0.001, A = 0.4, c = 0.1, A = 0.45, lu = 0.1, s(0) = 0, 0(0) = tt, and (a) Ai = 0.002, (b) Ai = 0.3 (c) and 
(d): The corresponding power spectrum, where the parameters in Fig. 6(c)are the same with Fig. 6(a) and the parameters in 
Fig. 6(d)are the same with Fig. 6(b). 
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